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Abstract 

This paper provides an alternate derivation of the equations used in the GRMHD 
code of De Villiers & Hawley using Stokes Theorem. This derivation places the pub- 
lished form of the discretized equations of GRMHD in a broader context, and suggests 
possible future directions for numerical work. 

1 Introduction 

The general relativistic magnetohydrodynamic (GRMHD) code of De Villiers & Hawley 
(2003; DH) has been used to study the accretion of matter onto black holes. It is a 
highly stable 3 + ID code that has been used to simulate in great detail the evolution 
of accretion disks in the spacetime of Kerr black holes. This code has established the 
ubiquity of energetic axial jets in such environments (De Villiers et al, 2003, 2005), 
and has shed light on the magnetic environment produced accretion disks (Hirose et 
al, 2004; De Villiers, 2006). 

The refinement of the code through testing and development is ongoing. One area 
of study that may help spur future improvements is the mathematical foundation of 
the GRMHD equations. There are many threads in the literature documenting the 
differential version of the equations. Misner, Thorne, & Wheeler (MTW; 1973) pro- 
vide a clear derivation of the equations of hydrodynamics. Anile (1989) discusses the 
extension to magnetohydrodynamics (MHD) in a special relativistic background. Haw- 
ley, Smarr, & Wilson (1984; HSW) treat hydrodynamics in a fixed background (black 
hole spacetime), while DH completes the derivation of the complete set of GRMHD 
equations. 



This paper presents a derivation of the GRMHD equations from the generahzed 
Stokes Theorem (MTW; see also Page &: Thorne, 1974), and was motivated by a 
derivation due to Clarke (1996) of the equations of non-relativistic MHD as discretized 
for the ZEUS code. Clarke provides an integral formulation of the equations of MHD 
using Stokes Theorem, and clearly shows the connection between the apparently sim- 
plistic finite-differencing techniques used in ZEUS (and also in the GRMHD code) to 
the notion of flux conservation embodied in Stokes Theorem. 



2 Conservation Laws of GRMHD 



The state of an ideal magnetized fluid on a curved background is described by a set of 
conservation laws. The equations in differential form are 

V^(/9C/^)=0 (continuity), (1) 
'SJ ^T^^ = (energy — momentum), (2) 
y^*ptiu^Q (induction), (3) 

where is the covariant derivative, p is the density, the 4- velocity, T^^ the energy- 
momentum tensor, F^^^ the electromagnetic field strength tensor, and its dual is 

"F^^^ = ^e^'^'^^F^^, (4) 

where e^^^^ = —{1/yJ—g) [fiuS^] is the contravariant form of the Levi-Civita tensor. 
The energy-momentum tensor is 



4 vr \ 4 / 



(5) 



where h = 1 + e + P/p is the specific enthalpy, with e the specific internal energy and 
P = pe (r — 1) the ideal gas pressure (F is the adiabatic exponent). 

Following Lichnerowicz (1967) and Anile (1989), define the magnetic induction and 
the electric field by projecting onto the four- velocity, 

B'' = *F^"'U^, (6) 
^ F'^'U^,. (7) 

The fiux-freezing condition, wherein the electric field in the fiuid rest frame is zero, 
implies F^^'^ Ui, = 0. Combine the definition of the magnetic induction ^ with (^J and 
flux freezing to obtain 

Ff,u = e^(,^,B"U'^, (8) 
where e^^sy = V~9 it^'^ ^l]- The orthogonality condition 

B''U^ = (9) 
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follows directly from ^ and the anti-symmetry of F^^'^ . 

Using these results, it is possible to rewrite the electromagnetic portion of the 
energy-momentum tensor to obtain 

Tf"" = (/9/i+||6f) U^U''+ ^^ + ^^ gf"" -bf'b'' (10) 

where ||6|p = b^b^ is the magnetic field intensity, and = B^^/{4tt) is the magnetic 
field four-vector. 

In addition, we make the identification 

B"^ = F^e,B' = Fr^,13''' = Fer, (11) 

where are the Constrained Transport (CT; Evans & Hawley, 1988) magnetic field 
variables. Also, 

(12) 



where are the CT EMFs. 



Ftr 


= V't'B'^ 


- B'^ 


= £." 


Fte 


= v"" B'f' 


- v'^ B"" 


= £' 


Ft4, 


= V^B' 


-VB^ 


= 8'^, 



3 Stokes Theorem 

The generalized form of Stokes Theorem (MTW) relates the divergence of a vector 
field (^'^) in a 4-volume (V) to the flux of this vector field through the oriented surface 
(9V) bounding this volume: 

/ V^,A^d^^= I A^cfj:^ (13) 

Jv JdV 

where d'^Q is the volume element of V and d^S^ is the four- vector outward normal 
to the surface dV. Since the GRMHD code uses the Boyer-Lindquist coordinates for 
the Kerr metric, the time and space coordinates (t, r, 6, (p) will be used from hereon. 
In Boyer-Lindquist coordinates, the volume element is d^VL = —g dt dr d6 dcj) where 
= a ^/7, a = (—5**) ^^"^ is the lapse function, and ^ is the determinant of the 
spatial 3-metric. 

Expression (|13j) is completely general; in anticipation of its application to the 
GRMHD code, take the four-volume V to be an infinitessimal region whose geome- 
try matches that of a discrete 3 + \D zone in the computational grid. This four- volume 
has eight well defined boundary planes, corresponding to hypersurfaces where one the 
the Boyer-Linquist coordinates is held constant. Figure ^ presents a three-dimensional 



3 



sketch of the bounding faces, and only the bounding faces in the radial direction are 
labelled. The surface denoted ro represents the hypersurface of constant r = ro, the 
"lower" bounding hypersurface of V in the radial direction. The surface denoted r^ + dr 
represents the hypersurface of constant r = ro + dr, the "upper" bounding hypersurface 
of V in the radial direction. Similarly, the hypersurfaces of constant polar angle are 
denoted Oq and 9o + d9; and hypersurfaces of constant azimuthal angle, 0o ^-nd (pQ + dcj). 
The generalized form of Stokes Theorem also deals with the time coordinate on an 
equal footing, so there are also surfaces of constant time, to and + dt-, bounding V. 




Figure 1: Sketch of hypersurfaces bounding the four-dimensional volume V. 



Figure ^ also introduces some additional terminology for the staggered mesh used 
in the GRMHD code: scalar quantities (e.g. p) reside on zone centres, whereas vector 
quantities (e.g. B'^) reside on boundary surfaces (zone faces). In addition, the discrete 
grid associates the face-centred quantities to the "lower" bounding hypersurface (e.g. 

associated with the zone of Figure ^ lies on the centre of the ro hypersurface, while 
the B^ on the ro + dr hypersurface is associated with the next zone in the radial 
direction) . 



4 Discretizing the Integral Equations 



The discretization of the integral form of the equations of GRMHD can be illustrated 
by considering a conservation law for some arbitrary vector field, V^^^ = 0, which 
can be expressed using Stokes Theorem (fT^ as an integral over dV: 

A^'d^T.^ = Q. (14) 

dV 

Using the geometry of volume V introduced in the previous section, the conservation 
law reads. 







/ A^'d^^^, = j A'' dt dO d<l) - j A'' dt dO d(t> (15) 

J dV J ro+dr J ro 
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A" yf^dtdrd(j)- I A" ^f^dtdrd<i) 

eo+de Jeo 

+ I A'^ ^/^dtdrde - I A'^^/^dtdrde 

J dn+dd> J (ho 



+ / A'V^drdOd^- / A'V^drdOd^ 

J to+dt J to 

where the subscript on each integral indicates the specific hypersurface (see Figure 
and each integral symbol in understood to represent a definite triple integral over each 
bounding hypersurface. 

The presence of the time coordinate in the integrals over surfaces of constant spatial 
coordinate serves as a pointed reminder that proper centering of physical quantities 
must be done not only in space but in time as well. If for no other reason, the natural 
emergence of this important connection between all four coordinates validates the use 
of the integral formulation as the starting point for numerical discretization. 

In going to a discrete grid, letting dx^ — > Ax^, (|15|) becomes an approximate 
relation, and there are many possible ways to proceed with the discretization. To see 
how this comes about, consider the fiux of A^ through the surface of constant tq, 

f f(f)o+M fOo+Ae fto+At ~\ 

A'' ^/^dtded(|)= I / / A'\/^dtd9d(l)\ (16) 

ro W</>0 Jeo J to J ro 

In subsequent expressions, the braces { } will denote a given hypersurface; this nota- 
tion is similar to that used by Page & Thorne (1974). Neglecting the integral over time 
for a moment, the simplest approximation to (|16j) is the product of the mean radial 
component of A on the surface multiplied by the area of the surface, 

<f>o+A<l> reo+Ae ~\ f r<i>o+A(j> rOo+Ae 

/ A'-^dedci^i -{A'lrol / V^dOd 

<f>0 Jeo J ro yJ't>0 J Go J ro 

(17) 

where {A'^)ro represents the mean value of A^ on the hypersurface rg. On a discrete 
grid, the fiux could be approximated by taking quantities on the zone face as repre- 
sentative of this mean so that 



C r4>o+A4> i-eo+Ae 

/ / A'V^ded^\ ^{A^V^AOAcI)}^^, 

kJ4>o J do J ro 



(18) 



where A^ and \/—g on the right-hand side are understood to be evaluated at the centre 
of the ro zone face. 

However, equation (|16j) also includes the integral over time which must be dealt 
with in a consistent manner. In a numerical implementation, the solution at time level 
to is known and the solution at time level to + is what is sought, so (|16|) cannot 
be evaluated exclusively with data at the known time level to without introducing 
potentially serious numerical errors. For consistency, (|TH|) should be evaluated in such 
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a way as to involve the unknown time level. For example, a complete discretization of 
((TB|) could be written as 



C /-eo+^f) i-tQ+m 1 1 r r i i 

{ / / A''^dtd9d^\ «-| jr{n) j^jr{n+l) ^AtA^A^l , 

U0O J do J to J,.(, 2 U J Jro 

(19) 

where solutions at different time levels are distinguished by introducing superscripts; 
^'■{") denotes the value of at the current time t^, and ^ at the next time 
level, to + At. 

Similar expressions can be derived for the other hypersurfaces of constant spatial 
coordinate. This leaves the two hypersurfaces of constant time. Take the surface of 
constant to, the flux through this surface is 



r C i-<f>o+^4> r0o+Ae pro+Ar ^ 

/ A^ dr dO del) = < / / A^^/^drded^\ 

J to \.J4>oJdoJro J 



(20) 



to 



which can be approximated as 

r>(/>o+A0 i-do+Ae j-ro+Ar 
Jro 



A*- y/^drdOt 



A 



tin) 



-gArAOAcp. (21) 



to 



Note that for the hypersurfaces of constant time, the brace notation is not used, and 
the time-components of vectors are taken to be zone-centred quantities. 

Combining all these terms, and, in keeping with the GRMHD code, taking At to 
be constant over the entire grid, at all times, ()15() can be approximated as 







+ 



^t (n+l) _ (n) 

— + 7^ 



At 



L L J J ro+Ar I L J J 



Ar 



}en+Ae { 



do 



AO 



; J ) 4IO + A4, I L J ) 410 

A(f> 



(22) 



Of course, this is but one possible discretization of ()15|). The advantage of the 
integral formulation is that it facilitates writing flux-conserving numerical routines 
which achieve proper centering not only in space but in time as well. When working 
with a differential formulation of the conservation laws, such choices may not be as 
obvious. 
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5 Equations of GRMHD in Integral Form 



The numerical solutions to the equations of GRMHD follow the time-evolution of a 
fluid on a 3D spatial grid from some prescribed initial state at some Boyer-Lindquist 
time ti. The numerical solver advances the solution in small, uniform increments At 
and produces a time-ordered sequence of states for the fluid. 

The integral form of the GRMHD equations is 



As in the previous sections, these equations will be discretized by taking the volume V 
to represent a zone on the computational grid. 

The previous section showed how an implicit scheme could be obtained from Stokes 
Theorem. Instead of using such a scheme, the GRMHD code approximates the tempo- 
ral integrals by building an approximate, upwinded solution at an intermediate time 
level; this approach gives rise to an explicit scheme that is computationally simpler. 
An additional motivation is that the conservation laws involve a mix of scalar and vec- 
torial code variables which must be properly centred. For example, to achieve proper 
centring in time and space in the case of the equation of continuity, the product of 
density p (zone-centred) and four-velocity (face-centred) is evaluated by upwinding 
the density, i.e. by interpolating density to the appropriate zone face and advancing it 
in time to the intermediate time level, n + 1/2. This estimate, at the half-step, is then 
used in the continuity equation. Denote the upwinded quantity using an overbar; for 
instance, radially upwinded density would be shown as p*" C/^. 

The details of the upwinding technique are not relevant to the present discussion, 
as they simply represent an explicit means of achieving proper time-centring of the 
integrals over the time coordinate in (|23|) . A description of the technique as it applies 
to the GRMHD code can be found in DH and in HSW. 




(23) 
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5.1 Continuity 



The integral form of the continuity equation is 

JdV 

Using the geometry of volume V introduced earlier, the equation reads, 
= / {pU'')^/^dtded(^- I {pU'')^/^dtd9d(j) 

J m+dr J m 



(24) 



(25) 



Oo+de 



{pU^)^f^dtdri 



{pU'^)^/^dtdri 



do 



{pU't')^/^dtdrdd- / {pU'f')^^dtdrde 
(f>o+d(f> J 4>o 

+ I {pU^)y/^drded(t)- I {pU^)./^drd9d(p. 

Jto+dt J to 

Using the upwinding notation we have, for example, 

[ (p U'') ^ dt de # « {r At AO A0}^^ . 



(26) 



Using the definitions of transport velocity, = W/a, and auxiliary density, D 
pW, this result is rewritten 



[ {pU'')^dtd9d<p^{D''V'^AtA9A(j>}^^. 
The integrals through the hypersurfaces of constant t are approximated as 



(27) 



/ {pU^) dr dO d(l) - I {pU^)./^drde, 

J to+dt J to 



Combining these results and simplifying yields 



^ArAOAcj) 
(28) 



At 



1 

7f 



ro+A r 



ro 



Ar 



(29) 



+ 



+ 



Oo 



A6 



which is the discretized, upwinded form of equation (24) in DH as implemented in the 
GRMHD code, recovered very simply here by the integral formulation of the equation 
of continuity. Note how this approach "breaks" the implicit coupling at the unknown 
time level on the left- and right-hand sides shown in (|22() . and achieves proper centering 
through an explicit method. 
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5.2 Induction 



The integral form of the induction equation is 

*F^"' d^S^ = 0. 

dV 



(30) 



Using the geometry of volume V introduced earlier, the four components of the induc- 
tion equation are, 







dV 



+ 



ro-\-dr 



do+dO 



-g dt dO ( 



^F^^^/^dtdn 



+ 



I *F^'' y/^dtdrde 
*F*''^/^drde d^- 



^F''" ^/^dtdedcj) (31) 
^F^" ,/^dtdrd(j) 
*F^'' ./^dtdrdO 



00 



to+dt 



* F'^" yf^ dr dO ^ 



to 



Though the previous examples did not deal with tensor quantities, the results are 
readily extended. Taking the time component, v = t, the integrals over the surfaces of 
constant t are identically zero (antisymmetry of F^"), leaving 







dV 



* rjrt 



F"^/^dtde, 



ro+dr 



'F 



gdtd9d(j) (32) 
g dt dr dcf) 

+ I *F't'^,/^dtdrde- I *F'i'^ ./^dtdrdO. 



+ 



*F^^ ^f^dtdrdcj) 



Oo+de 



do 



Using ^ and substituting the expressions for the CT variables, ((TT|) and (fT^ . back 
into (ISH) yields: 







B'^dtdOi 



' ro+dr 

+ / B^dtdrt 



B'' dt de i 



(33) 



B'^ dtdr. 



do 



+ / B'^dtdrde- / B'f'dtdrde. 

'(t)o+d4> J 4)0 

To discretize this result, note that the vanishing of the integrals over hypersurfaces of 
constant t indicates that this relation must be satisfied at all times; this component 
of the induction equation is a constraint on the CT magnetic field. Equation (|33jl 
simplifies to: 







mr,^Ar-mr, , {^'},„+A. " {^'l.o , W 



Ar 



+ 



(f>0+^<l> 



{B^} 



<t>0 



(34) 
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which is the famihar "div B" constraint on the magnetic field, and equivalent to equa- 
tion (47) of DH. 



The spatial components of the induction equation describe the time-evolution of the 
CT magnetic field. To see how this comes about, take the radial component, v = r, 
of the induction equation. The integrals over the surfaces of constant r are identically 
zero, leaving 







'F^'\f^dtdr d(l) 



(35) 



+ / *F'f"' ./^dtdrde - I *F'^' 



-g dt dr dO 



+ 



* F^"^ ./^ dr dO ^ 



to+dt 



*F^'' ^/^dr dOi 



to 



Using Q), ((TT|) . and (fT^ . equation can be discretized and rearranged to read 



At AO A(j) 

Similar rearrangements can be made for the 9 and (p components of (|31jl 



(36) 







At 



^"^1 -1^*1 iri -in 

Ar A(j) 



_ ^^{n) ^£ ir-p+Ar l*^ irp _ Igp+Ag Igp 

At Ar A6I 



(37) 



(38) 



These expressions are formally equivalent to the CT equations of DH (cf. equa- 
tion (34) in DH). However, as with the continuity equation, the issue of proper time- 
centering, this time of the EMFs, is important. This is flagged in the above expressions 
by using an overbar, £, on the EMFs. 

Evans & Hawley (1988) discussed the issue of time-centering of the EMFs, noting 
that there seemed to be some flexibility to the prescription of the EMFs. DH provided 
a derivation of a numerical formulation similar to the MOCCT approach of Hawley 
&: Stone (1995) to provide time-centred expressions which were stable under a variety 
of numerical tests. Perhaps the main advantage of revisiting CT from the point of 
view of Stokes Theorem is that this formulation makes it more obvious that proper 
centering, both in space and time, is indeed required. A prescription for EMFs that 
is compatible with the integral formulation has the EMFs edge-centred on the spatial 
grid and centred on the half-step on the temporal grid. The approximation using 
incompressible MHD documented in DH is but one possible realization of this. Of 
course, implicit formulations of (|^T|l are also possible, but are likely to be difficult to 
implement. 
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5.3 Energy-Momentum Conservation 



The integral form of the energy- momentum conservation equations, 

T^'" d^S. = (39) 



IdV 

Using the geometry of volume V introduced earlier, the four components of the energy- 
momentum equation are, 

0= /" T^^d^^f, = I ^dtded<t)- I T''" ^dtded<t) (40) 

J dV J ro+dr J ro 



T"" ^/^dtdrd4>- / T"" ^/^dtdrdct) 
+ / T^" ^/^dtdrde - I T'^" ./^dtdrde 



+ / dr dO d(t) - / dr de d(p. 

Jto+dt J to 

The discretization of this set of equations proceeds as described above. The only 
issue is the relation between the primitive code variables and T^'^ , and the extraction 
of the primitive variables. 

In the GRMHD code, the energy-momentum conservation law is treated differently, 
and the resulting equations provide a so-called non-conservative treatment of energy. 
The conservation law V^T'^'^ = is applied to a fluid with four-velocity W^. A very 
convenient reformulation of the conservation law arises when one "reorients" the equa- 
tions into a component parallel to the local four-velocity and components orthogonal 
to it (MTW). Projection onto the four- velocity yields a scalar relation, 

U,V^T^"' = 0, (41) 

which, as shown in HSW and DH, leads to the equation of internal energy conservation, 

dt{E) + ^di{^EV')+Pdt{W) + ^d^{^WV') =0. (42) 

To obtain the components of the conservation law orthogonal to the four-velocity, use 
the projection tensor P^^, = g^^ + Ufj,Uu, 

Pau T'^" = ga. T''" + T'^^ = T'^^ = 0, (43) 

where the final result was obtained by commuting the metric with the covariant deriva- 
tive and also using (|^T|) . As shown in DH, the three spatial components of this equation 
are 

= 5t {Sj - a hj 6*) (44) 
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The time-component is not used; instead, the GRMHD code uses the normahzation 
condition 



This paper has outhned the development of discretized equations of ideal General Rel- 
ativistic Magnetohydrodynamics (GRMHD) from the integral form of the conservation 
laws. Stokes Theorem transforms the conservation laws of GRMHD into simple flux 
integrals through the boundary of a closed four-volume, and many possible discretiza- 
tions are possible for these integrals. The particular choices used in the GRMHD code 
for the equation of continuity and the induction equation represent but one implemen- 
tation, arguably the simplest and most straightforward of these. Many other possible 
implementations are implied by equations H25() and H31|) . The GRMHD code uses an 
alternate discretization of the energy-momentum conservation equation; the integral 
formulation discussed here, (|4U|) . offers another possible discretization in keeping with 
the treatment of the continuity and induction equations. A note of caution is in order: 
the integral relations derived here, dealing as they do with coupled sets of non-linear 
equations, may not yield stable numerical schemes. Since only limited stability analysis 
is possible in non-linear differential equations, theoretical work must be complemented 
by careful development and testing to validate any particular numerical treatment. 

In addition, this discussion has implicitly assumed smooth flows; the issue of shock 
capturing has not been touched upon. The GRMHD code uses a relativistic treatment 
that is based on the original Von Neumann &i Richtmyer (1950) artificial viscosity 
treatment, a relativistic extension of this technique by May &: White (1967), and HSW. 
This remains the state of the art as of this writing. 

As a final note, the zone-based, flux-conserving formulation discussed here may also 
be of benefit in treating boundary conditions. The boundary zones, or ghost zones, 
of the computational grid are used to impose physical boundary conditions on the 
computational volume and also to provide approximate data used by the interpolation 
routines to provide correct time-centering of zones on the computational grid adjacent 
to the physical boundaries. The main interest in the integral formulation of boundary 
data lies at the inner radial boundary, where frame dragging effects are important, 
and also at the difficult "corner zones" which sit just outside the pole caps of the 
black hole. Definitive statements regarding the flux of electromagnetic energy near 
the event horizon hinge on reliable boundary treatments, and it is hoped that a more 
sophisticated approach, based on the integral formulation discussed here, may help. 
This remains an open area of research. 




(45) 



to close the set of equations. 



6 Summary and Discussion 
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